*This program runs a probit and calculates the associated marginal effects
*for stocks and the house, incorporating pension wealth
*It's used on the paper on wealth counterfactuals 
*for the verification for the Review of Economics and Statistics
*The program is run on Stata 11, using the version 9.2 option
*June 2011 - Dimitris Christelis

version 9.2
cap clear all
macro drop _all
clear
set mem 500m
set more off

*Paths
do "c:/foldp/foldp.do"
global projf "${foldp}/Wealth Counterfactuals/Verification"

*Tolerance
sca tolr=1e-6

*Weights
global wgt wgth

*ATTENTION: The base country (Germany in our case) has to be put first in the sequence
global eureg1 uk
global eureg2 England
*Country numbers: they have to follow exactly the same sequence as the country names
global eureg3 12   
local nas: word count $eureg1

global assetl1 rfin home /*ownb secdebt liab tdebt*/
global assetl2 hrfin hhome /*hownb hsecdebt hliab htdebt*/
local nar: word count $assetl1


*Attention: the asset loop has to be before the country loop
forvalues ll=1/`nar' {

   *Globals of asset
   global as: word `ll' of $assetl1
   global hhas: word `ll' of $assetl2

   forvalues kk=1/`nas' {

   *Globals of country
   global re: word `kk' of $eureg1
   global fnre: word `kk' of $eureg2
   global cn: word `kk' of $eureg3


   *Loop for inclusion of pension wealth   
   forvalues www=1/3 {   
    
   if `www'==1 { 
      global pw nopw
      global aw 0
   }
   if `www'==2 { 
      global pw pw0
      global aw hpw0
   }
   if `www'==3 { 
      global pw pw1
      global aw hpw1
   }

   *Logging
   cap log close
   log using "${projf}/Results/ELSA/penwealth/${as}/penwealth_${re}_${as}.log", replace


   global pdate "$S_DATE" 
   global ptime "$S_TIME"
   display "log file for ${as} in ${fnre} opened on ${pdate} at ${ptime}"

   *Opening HRS data to generate quartiles
   use "${projf}/Data/HRS/HRS_DG.dta", clear
   
   *Keeping only the pension wealth sample
   qui keep if pwsample==1

   *Putting pension wealth with a 1% real growth rate of pension benefits
   *equal to the one with 0% growth, so that it can be compared with the 
   *European countries
   qui gen double hpw1 = hpw0
    
   *Checking the sample weights
   assert abs(wgth-wgtach)<1e-6
   drop wgtach

   *Keeping only heads
   keep if head==1 & wgth>0 & wgth!=.

   *Total debt
   qui gen double htdebtv = hmortv+hhmlov+hliabv
   qui gen byte htdebto = (htdebtv>tolr)

   *Secured debt, except on 2nd house
   qui ge double hsecdebtv = hmortv + hhmlov /*+ hmort2v*/
   qui ge double hsecdebto = (hsecdebtv>tolr)

   *Creating net worth and income quartiles
   *Net worth, equal total minus the asset
   if "${as}"=="rfin" | "${as}"=="home" | "${as}"=="ownb" {
      qui gen double nw=hnetwv-${hhas}v+${aw}
   }
   *For debts, we add back instead of substracting
   if ! ("${as}"=="rfin" | "${as}"=="home" | "${as}"=="ownb") {
      qui gen double nw=hnetwv+${hhas}v+${aw}
   }
   
   *Income, equal to total minus capital income
   qui gen double inc=hgtincv-hcpincv

   *Creating wealth quartiles 
   qui su nw if head==1 [aw=wgth], detail
   sca p25_w_us = r(p25)
   sca p50_w_us = r(p50)
   sca p75_w_us = r(p75)

   *Creating income quartiles
   qui su inc if head==1 [aw=wgth], detail
   sca p25_i_us = r(p25)
   sca p50_i_us = r(p50)
   sca p75_i_us = r(p75)

   *Reading the SHARE data
   use "${projf}/Data/SHARE/SHARE_DG.dta", clear

   *Keeping only the country of interest
   qui keep if country==3 & ${wgt}>0 & ${wgt}!=.
   
   *Keeping only the pension wealth sample
   qui keep if pwsample==1

   *Keeping only heads
   qui keep if head==1

   *Total debt
   qui gen double htdebtv = hmortv+hliabv /* +hhmlov*/
   qui gen byte htdebto = (htdebtv>tolr)

   *Secured debt, except on 2nd house
   qui ge double hsecdebtv = hmortv /*+ hhmlov + hmort2v*/
   qui ge double hsecdebto = (hsecdebtv>tolr)

   *Creating net worth and income quartiles
   *Net worth, equal total minus the asset
   if "${as}"=="rfin" | "${as}"=="home" | "${as}"=="ownb" {
      qui gen double nw=hnetwv-${hhas}v+${aw}
   }
   *For debts, we add back instead of substracting
   if ! ("${as}"=="rfin" | "${as}"=="home" | "${as}"=="ownb") {
      qui gen double nw=hnetwv+${hhas}v+${aw}
   }
   
   *Income, equal to total minus capital income
   qui gen double inc=hgtincv-hcpincv

   qui su nw if head==1 [aw=wgth], detail
   sca p25_w_de = r(p25)
   sca p50_w_de = r(p50)
   sca p75_w_de = r(p75)

   qui su inc if head==1 [aw=wgth], detail
   sca p25_i_de = r(p25)
   sca p50_i_de = r(p50)
   sca p75_i_de = r(p75)

   *Reading the parameters, variable lists etc.
   do "${projf}/Programs/Parameters.do"

   *global nsim=2
   *global finsime=5*${nsim}

   use "${projf}/Data/ELSA/ELSA_DG.dta", clear

   *Keeping only the country of interest
   qui keep if country==${cn} & ${wgt}>0 & ${wgt}!=.

   *Keeping only heads
   qui keep if head==1

   *Keeping only the pension wealth sample
   qui keep if pwsample==1

   *Total debt
   qui gen double htdebtv = hmortv+hhmlov+hliabv
   qui gen byte htdebto = (htdebtv>tolr)

   *Secured debt, except on 2nd house
   qui ge double hsecdebtv = hmortv + hhmlov /*+ hmort2v*/
   qui ge double hsecdebto = (hsecdebtv>tolr)

   *Ownership dummy
   *qui gen byte ${hhas}o=(${hhas}v>tolr)
   
   *Rescaling age
   qui replace age=age/sc_age
   qui gen double agesq=age^2

   *Creating net worth and income quartiles
   *Net worth, equal total minus the asset
   qui gen double nw=hnetwv-${hhas}v+${aw}
   *Income, equal to total minus capital income
   qui gen double inc=hgtincv-hcpincv

   *Bad health dummy
   qui gen byte badhs=(srhealtha>=4)

   *College dummy
   qui gen byte college=(edu==3)
   
   *High school dummy
   qui gen byte highs=(edu==2)

   *Couple dummy
   cap drop couple
   qui gen byte couple=(mstatus==1)
   
   *Widow dummy
   qui gen byte widow=(mstatus==3)
   
   *Never married
   qui gen byte nevmar=(mstatus==4)
      
   *Creating wealth quartiles for financial and real assets
   qui gen double wqnw1_qus=(nw<=scalar(p25_w_us))
   qui gen double wqnw2_qus=(nw>scalar(p25_w_us) & nw<=scalar(p50_w_us))
   qui gen double wqnw3_qus=(nw>scalar(p50_w_us) & nw<=scalar(p75_w_us))
   qui gen double wqnw4_qus=(nw>scalar(p75_w_us))
   qui gen double wqnw1_qde=(nw<=scalar(p25_w_de))
   qui gen double wqnw2_qde=(nw>scalar(p25_w_de) & nw<=scalar(p50_w_de))
   qui gen double wqnw3_qde=(nw>scalar(p50_w_de) & nw<=scalar(p75_w_de))
   qui gen double wqnw4_qde=(nw>scalar(p75_w_de))

   qui gen double wqinc1_qus=(inc<=scalar(p25_i_us))
   qui gen double wqinc2_qus=(inc>scalar(p25_i_us) & inc<=scalar(p50_i_us))
   qui gen double wqinc3_qus=(inc>scalar(p50_i_us) & inc<=scalar(p75_i_us))
   qui gen double wqinc4_qus=(inc>scalar(p75_i_us))
   qui gen double wqinc1_qde=(inc<=scalar(p25_i_de))
   qui gen double wqinc2_qde=(inc>scalar(p25_i_de) & inc<=scalar(p50_i_de))
   qui gen double wqinc3_qde=(inc>scalar(p50_i_de) & inc<=scalar(p75_i_de))
   qui gen double wqinc4_qde=(inc>scalar(p75_i_de))

   *Saving the dataset
   qui save "${projf}/Results/ELSA/penwealth/${as}/temp.dta", replace

   *Running 2 regressions, using both definitions of quartiles
   foreach cb in qus /*qde*/ {

      *Bootstrap program

      *Creating a matrix with number of columsn equal to 1 + # bootstrap runs, storing
      *the main sample estimates in the first line, together with the proportion of
      *hhds owning the asset
      qui su ${hhas}o if head==1 [aw=$wgt], meanonly
      sca ${as}p=r(mean)
      mat ${as}coeffall=${as}p
      local nch=$nc+1
      mat colnames ${as}coeffall = _b`nch'   
      qui save "${projf}/Results/ELSA/penwealth/${as}/temp2.dta", replace
  
      *Subsequent runs, beginning of the bootstrap loop
      forvalues i=1/$finsim {

         if `i'==1 {
            local sr = 1
         }
         *Augmenting the index that tracks valid bootstrap runs by one
         if `i'>1 {
            local sr = `sr' + 1
         }
           
         use "${projf}/Results/ELSA/penwealth/${as}/temp2.dta", clear
         *Fixing the seed before each resampling
         local sd=1234+30*`i'+`kk'
         set seed `sd'
         bsample
               
            qui su ${hhas}o if head==1  [aw=$wgt], meanonly
            sca ${as}p=r(mean)
            mat ${as}coeff`sr'=${as}p
            mat ${as}coeffall=(${as}coeffall \ ${as}coeff`sr')

            forvalues lll=1/10 {
               if `sr'==$nsim*(`lll'/10) {
                  noisily di in yellow `sr' in yellow "*" _continue 
               }
            }

            mat drop ${as}coeff`sr'

         
               
         *If the number of valid bootstrap runs is equal to the predetermined
         *one then the loop over bootstrap runs stops
         if `sr'==$nsim {
            di in yellow "End of ${nsim} valid boostrap runs, total runs equal to `i'"
            continue, break
         }
         
      }  /* End of loop over bootstrap runs */

      *Saving the matrix of sample proportions
      drop _all
      svmat double ${as}coeffall
      qui save "${projf}/Results/ELSA/penwealth/${as}/boot_predpr_${re}_${as}_${pw}_`cb'.dta", replace

        
      ***********************************************************
      *Calculation of marginal effects and associated magnitudes*
      ***********************************************************

      global nc: word count ${basiccovs_`cb'} _cons

      *Reading the dataset
      drop _all
      use "${projf}/Results/ELSA/penwealth/${as}/temp.dta", clear
      keep if pwsample==1
      
      sca nobs_`cb'=_N

      *Incrementing the continuous variables by the appropriate magnitudes
      *to compute marginal effects and elasticities
      *Variables that are not used in elasticity form
      foreach y in $conlist_meff {
         qui gen double `y'_pl=`y'+(ch_`y'/sc_`y')
      }   

      *Merging with the sample proportions
      merge using "${projf}/Results/ELSA/penwealth/${as}/boot_predpr_${re}_${as}_${pw}_`cb'.dta"
      drop _merge
      
      *Renaming the coefficient of the constant 
      *set trace on
      rename ${as}coeffall pr_${re}_${as}_${pw}_`cb'_own

   
      *Merging with the draws from the US
      if "`cb'"=="qus" {
         merge using "${projf}/Results/HRS/Country/penwealth/${as}/boot_drawcoeff_usall_${as}_${pw}.dta"   
      	 drop _merge   

      	 *Putting the values of the drawn coefficients equal to the dummy value from missing
      	 *in rows lower than the number of bootstraps, so that there are no missing values
      	 local nch = $nc + 1
      	 forvalues h=1/`nch' {   
      	    qui replace _b`h'=$dumval if _n>${nsim}+1  
      	 }      
   
      	 *Renaming the coefficients  
      	 local ncl = $nc - 1
      	 forvalues h=1/`ncl' {   
      	    local g: word `h' of ${basiccovs_qus}
      	    rename _b`h' b_usall_`g'   
      	 }   

      	 *Renaming the coefficient of the constant 
      	 *set trace on
      	 rename _b${nc} b_usall_cons
   
      	 local nch = $nc + 1
      	 rename _b`nch' pr_usall_${as}_${pw}_own
      }
      
      *Merging with the coefficients and the probabilities from Germany
      if "`cb'"=="qde" {
         if "${re}"~="de" {

      	    *Coefficients
      	    merge using "${projf}/Results/SHARE/de/${as}/boot_drawcoeff_de_${as}_${pw}_qde.dta"   
      	    drop _merge   

      	    *Putting the values of the drawn coefficients equal to the dummy value from missing
      	    *in rows lower than the number of bootstraps, so that there are no missing values
      	    local nch = $nc + 1
      	    forvalues h=1/`nch' {   
      	       qui replace _b`h'=$dumval if _n>${nsim}+1  
      	    }      
   
      	    *Renaming the coefficients  
      	    local ncl = $nc - 1
      	    forvalues h=1/`ncl' {   
      	       local g: word `h' of ${basiccovs_qde}
      	       rename _b`h' b_de_`g'_qde   
      	    }   

      	    *Renaming the coefficient of the constant 
      	    *set trace on
      	    rename _b${nc} b_de_cons_qde
   
      	    local nch = $nc + 1
      	    rename _b`nch' pr_de_${as}_${pw}_qde_own
      	 }         
      }
      
      *Loop over bootstrap coefficients
      forvalues i=1/$finsim { 

         *Defining an index of valid bootstrap runs that will determine whether the 
         *loop finishes or not
         if `i'==1 {
            local sr = 1
         }
         *Augmenting the index that tracks valid bootstrap runs by one
         if `i'>1 {
            local sr = `sr' + 1
         }
         *If there were any condition on the coefficients that was violated by the current
         *draw then the index of valid bootstrap runs would be decreased by one at this point
         *in the program. For the pooled probit there is no such condition and thus the
         *index of valid bootstrap runs (sr) is equal to the on of total bootstrap runs (i)
   
      
         *Creating a counterfactual prediction using the coefficients from the US
         if "`cb'"=="qus" {
            qui gen double xb0_usall = b_usall_cons[`i']
            foreach y in $basiccovs_qus {
                qui replace xb0_usall = xb0_usall + b_usall_`y'[`i']*`y' 
            }
	 }

         *Creating a counterfactual prediction using the coefficients from Germany
         if "`cb'"=="qde" {
            if "${re}"~="de" {         
               qui gen double xb0_de_qde = b_de_cons_qde[`i']
               foreach y in $basiccovs_qde {
                  qui replace xb0_de_qde = xb0_de_qde + b_de_`y'_qde[`i']*`y' 
               }
            }            
         }
      
         *******************************************************************
         *Calculating the average of probabilities across individual obs
         if "`cb'"=="qus" {   
            qui gen double dp_usall=normal(xb0_usall)
            qui su dp_usall [aw=$wgt], meanonly
            sca r0m=r(mean)
         }
         
         *Probability using Germany coefficients
         if "`cb'"=="qde" {
            if "${re}"~="de" {
               qui gen double dp_de_qde=normal(xb0_de_qde)
               qui su dp_de_qde [aw=$wgt], meanonly
               sca r2m_qde=r(mean)
            }
         }
         if `i'==1 {
             
            if "`cb'"=="qus" {
               *Counterfactual probability using US coefficients
               qui gen double pr_${re}_${as}_${pw}_qus_cus=r0m
             
               *Difference between US ownership proportion and
               *own ownership proportion
               qui gen double dpus_${re}_${as}_${pw}_qus_tot=pr_usall_${as}_${pw}_own - ///
                  pr_${re}_${as}_${pw}_qus_own
             
               *Difference between US and own probability calculated using
               *US coefficients (characteristics effect)
               qui gen double dpus_${re}_${as}_${pw}_qus_char=pr_usall_${as}_${pw}_own[1] - ///   
                  pr_${re}_${as}_${pw}_qus_cus[1]
             
               *Difference between own probability calculated using
               *US coefficients and own probability (coefficient effect)
               qui gen double dpus_${re}_${as}_${pw}_qus_coeff=pr_${re}_${as}_${pw}_qus_cus[1] - ///   
                  pr_${re}_${as}_${pw}_qus_own[1]
            } 

            if "`cb'"=="qde" {
               if "${re}"~="de" {
            
             	  *Counterfactual probability using German coefficients
             	  qui gen double pr_${re}_${as}_${pw}_qde_cde=r2m_qde
             
             	  *Difference between mean German predicted probability and
             	  *mean own predicted probability
             	  qui gen double dpde_${re}_${as}_${pw}_qde_tot=pr_de_${as}_${pw}_qde_own - ///
             	     pr_${re}_${as}_${pw}_qde_own
             
             	  *Difference between German and own probability calculated using
             	  *German coefficients (characteristics effect)
             	  qui gen double dpde_${re}_${as}_${pw}_qde_char = pr_de_${as}_${pw}_qde_own[1] - ///
             	     pr_${re}_${as}_${pw}_qde_cde[1]  
             
             	  *Difference between own probability calculated using
             	  *German coefficients and own probability (coefficient effect)
             	  qui gen double dpde_${re}_${as}_${pw}_qde_coeff = pr_${re}_${as}_${pw}_qde_cde[1] - ///
             	    pr_${re}_${as}_${pw}_qde_own[1]
               }
            }
         }
         
         *Putting the probability at the line indexed by
         *the valid bootstrap run (sr), and not the general
         *run (i)
         
         if `i'>1 {
            if "`cb'"=="qus" {
               
               *Counterfactual probability using US coefficients   
               qui replace pr_${re}_${as}_${pw}_qus_cus=r0m in `sr'   
               assert pr_${re}_${as}_${pw}_qus_cus<$dumval if _n==`sr'   

               *Difference between US and own probability calculated using   
               *US coefficients (characteristics effect)   
               qui replace dpus_${re}_${as}_${pw}_qus_char=pr_usall_${as}_${pw}_own[`sr'] - ///   
                  pr_${re}_${as}_${pw}_qus_cus[`sr'] in `sr'   

               *Difference between own probability calculated using   
               *US coefficients and own probability (coefficient effect)   
               qui replace dpus_${re}_${as}_${pw}_qus_coeff=pr_${re}_${as}_${pw}_qus_cus[`sr'] - ///   
                  pr_${re}_${as}_${pw}_qus_own[`sr'] in `sr'   
            }

            if "`cb'"=="qde" {
               if "${re}"~="de" {
                  *Counterfactual probability using German coefficients
               	  qui replace pr_${re}_${as}_${pw}_qde_cde=r2m_qde in `sr'
               	  assert pr_${re}_${as}_${pw}_qde_cde<$dumval if _n==`sr'
         
               	  *Difference between German and own probability calculated using
               	  *German coefficients (characteristics effect)
               	  qui replace dpde_${re}_${as}_${pw}_qde_char=pr_de_${as}_${pw}_qde_own[`sr'] - ///
               	     pr_${re}_${as}_${pw}_qde_cde[`sr'] in `sr'

               	  *Difference between own probability calculated using
               	  *German coefficients and own probability (coefficient effect)
               	  qui replace dpde_${re}_${as}_${pw}_qde_coeff=pr_${re}_${as}_${pw}_qde_cde[`sr'] - ///
               	     pr_${re}_${as}_${pw}_qde_own[`sr'] in `sr'
               } 
            }
         }
         
         if "`cb'"=="qus" {
            drop dp_usall xb0_usall
         }
         
         if "`cb'"=="qde" {
            if "${re}"~="de" {
               drop dp_de_qde xb0_de_qde
            }   
         }   
      
         forvalues lll=1/10 {
            if `sr'==$nsim*(`lll'/10) {
               di in yellow "Iteration `sr', ${as}"
            }
         }

         *If the number of valid bootstrap runs is equal to the predetermined
         *one then the loop over bootstrap runs stops
         local ncpp=$nsim+1
         if `sr'==`ncpp' {
            di in yellow "End of ${nsim} valid bootstrap runs, total runs equal to `i'"
            continue, break
         }
      
      }  /* End of loop over bootstrap runs */

      *Dropping obs beyond the line of the last valid bootstrap run
      drop if _n>$nsim+1
   
      *Global of marginal effects, putting first the probabilities of interest
      if "`cb'"=="qus" {
         local me pr_${re}_${as}_${pw}_qus_own pr_${re}_${as}_${pw}_qus_cus dpus_${re}_${as}_${pw}_qus_tot ///
            dpus_${re}_${as}_${pw}_qus_coeff dpus_${re}_${as}_${pw}_qus_char
      }
      if "`cb'"=="qde" {
         if "${re}"=="de" {
            local me pr_${re}_${as}_${pw}_qde_own
         } 
         if "${re}"~="de" {
            local me pr_${re}_${as}_${pw}_qde_own pr_${re}_${as}_${pw}_qde_cde dpde_${re}_${as}_${pw}_qde_tot ///
               dpde_${re}_${as}_${pw}_qde_coeff dpde_${re}_${as}_${pw}_qde_char
         } 
      }
      

      *Keeping only the marginal effects and percentage variation
      keep country `me' `perc' 

      *Calculating the variance of the generated average probabilities and marginal effects
      *One has to start from the 2nd line, since the first line gives the magnitudes
      *evaluated at the estimated coefficients
      local ncpp=${nsim}+1
      foreach vr in `me' `perc' {
         qui su `vr' in 2/`ncpp' 
         qui gen double se_`vr'=sqrt(r(Var))
      }

      *Saving the dataset with all the probabilities and marginal effects
      qui save "${projf}/Results/ELSA/penwealth/${as}/boot_alldraweff_${re}_${as}_${pw}_`cb'.dta", replace

      *Saving the dataset with only the mean magnitudes and their standard errors
      qui keep in 1
      qui save "${projf}/Results/ELSA/penwealth/${as}/boot_meaneff_${re}_${as}_${pw}_`cb'.dta", replace
 
      *Saving matrices of marginal effects and percentage variations as a text file

      foreach rs in $prfx  {
         *Number of elements
         local n_`rs': word count ``rs''
         matrix olout=J(`n_`rs'',6,.)

         forvalues k=1/`n_`rs'' {
            local h: word `k' of ``rs''
            matrix olout[`k',1] = `h'[1]
            matrix olout[`k',2] = se_`h'[1]
            matrix olout[`k',3] = olout[`k',1]/olout[`k',2]
            matrix olout[`k',4] = 2*ttail(nobs_`cb' - ${nc},abs(olout[`k',3]))
            matrix olout[`k',5] = olout[`k',1] - invttail(nobs_`cb' - ${nc},sigl/2)*olout[`k',2]
            matrix olout[`k',6] = olout[`k',1] + invttail(nobs_`cb' - ${nc},sigl/2)*olout[`k',2]
         } 
    
         *Labels of matrix of magnitudes
         *Column labels
         local colmat `rs' se tstat pvalue ci_low ci_high
         matrix colnames olout = `colmat' 
         *Row labels
         matrix rownames olout = ``rs''
         *Displaying and saving the matrix of results
         mat list olout,format(%15.5f)
         mat2txt, mat(olout) ///
            saving("${projf}/Results/ELSA/penwealth/${as}/boot_mefftab_${re}_${as}_${pw}_`cb'.txt") replace format(%15.5f) ///
            title("Meff for ${as} in ${fnre}, using `cb' quartiles")

      }  /* End of loop over magnitudes */        
      erase "${projf}/Results/ELSA/penwealth/${as}/temp2.dta" 
   } /* End of loop over quartile specifications */
  
   erase "${projf}/Results/ELSA/penwealth/${as}/temp.dta"         

   }  /* End of loop over the kinds of pension wealth */

   }  /*End of loop over countries */    
}  /*End of loop over assets */    

display "log file for ${as} in ${fnre} opened on ${pdate} at ${ptime}"
display "log file for ${as} in ${fnre} closed on $S_DATE at $S_TIME"
cap log close
   
   
  


